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ABSTRACT 

We present a simplified solution to orbit determination of a binary system 
from astrometric observations. An exact solution was found by Asada, Akasaka 
and Kasai by assuming no observational errors. We extend the solution consider- 
ing observational data. The generalized solution is expressed in terms of elemen- 
tary functions, and therefore requires neither iterative nor numerical methods. 

Subject headings: celestial mechanics — astrometry — methods: analytical 
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1. Introduction 

In the seventeenth century, Kepler discovered the laws of the motion of celestial objects. 
Since then, a fundamental problem has been how to determine the orbital elements (and the 
mass) of a binary from observational data of its positions projected onto a celestial sphere. 
In fact, astrometric observations play an important role in astronomy through determining 
a mass of various unseen celestial objects currently such as a massive black hole (Miyoshi 
et al. 1995), an extra-solar planet (Benedict et al. 2002) and two new satellites of Pluto 
(Weaver et al. 2006). 

The orbit determination of resolved double stars was solved first by Savary in 1827, 
secondly by Encke 1832, thirdly by Herschel 1833 and by many authors including Kowalsky, 
Thiele and Innes (Aitken 1964 for a review on earlier works; for the state-of-the-art 
techniques, e.g, Eichhorn and Xu 1990, Catovic and Olevic 1992, Olevic and Cvetkovic 
2004). Here, resolved double stars are a system of two stars both of which can be seen. 
The relative vector from the primary star to the secondary is in an elliptic motion with a 
focus at the primary. This relative vector is observable because the two stars are seen. On 
the other hand, an astrometric binary is a system of two objects where one object can be 
seen but the other cannot like a black hole or a very dim star. In this case, it is impossible 
to directly measure the relative vector connecting the two objects, because one end of the 
separation of the binary, namely the secondary, cannot be seen. The measures are made in 
the position of the primary with respect to unrelated reference objects (e.g. quasars )whose 
proper motion is either negligible or known. 

There are two major problems in determination of a binary orbit: (1) observational 
errors and (2) the conditional equations which connect the observable quantities with the 
orbital parameters. The conditional equations are not only non-linearly coupled but also 
transcendental because of Kepler's equation (e.g., Goldstein 1980, Danby 1988, Roy 1988, 
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Murray and Dermott 1999, Beutler 2004). 

As a method to determine the orbital elements of a binary, an analytic solution in 
an exphcit form has been found by Asada, Akasaka and Kasai (2004, henceforth AAK) 
by assuming no observational errors. This solution is given in a closed form by requiring 
neither iterative nor numerical methods. The purpose of this brief article is to extend the 
AAK solution considering observational data. We will summarize the AAK solution in § 2, 
next extend the solution considering observational data in § 3, and make a computational 
test of the extended solution in 5 4. 



Here we briefly summarize our notation and formulation for ideal cases without any 
observational error. We consider only the Keplerian motion of a star around the common 
center of mass of a binary system by neglecting motions of the observer and the common 
center in our galaxy. Let {x, y) denote the Cartesian coordinates on a celestial sphere. A 
general form of an ellipse on the celestial sphere is 



which is specified by five parameters {a, l3,j,S,s) since the center, the major/minor axes 
and the orientation of the ellipse are arbitrary. Here, we should note that the origin of the 
coordinates is arbitrary in this paper, while it is taken at the position of the primary star 
in the case of resolved double stars (for instance, Aitken 1967). 

At least five observations enable us to determine the parameters by using the least 
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AAK solution 



ax^ + (3f + 2'yxy + 26 x + 2ey ^ 1, 
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square method as 
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where the location of the star on the celestial sphere at the time of tj for j = 1, ■ • ■ , n is 
denoted by Pj = {xj,yj). The inverse of a matrix in the L. H. S. exists uniquely. In some 
case such as a large observational error, however, the determinant of the matrix can become 
nearly zero. Then, it would be more difficult to estimate the inverse matrix. 

The time interval between observations is denoted by T{k,j) = tfc — tj for tk > tj. We 
choose the Cartesian coordinates (x, y) so that the observed ellipse given by Eq. ([T]) can be 
reexpressed in the standard form as 



0? IP' ' 
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where a >b. The ellipticity e is y 1 — IP/a^. There exists such a coordinate transformation 
as a combination of a translation and a rotation, which is expressed as 
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where p and a are given by 





(5) 



and VL is determined by 



tan 2Vt 
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Then, a and b are given by 



where / = ap^ + ^a^ + 2jpa, A = acos^fl + /?sin^O + 7 sin 20, and B = asin^O + 
/9cos^Q — 7sin2Q. 

It has been recently shown (Asada et al. 2004) that the orbital elements can be 
expressed explicitly as elementary functions of the locations of four observed points and 
their time intervals if astrometric measurements are done without any error. The key 
thing is that even after a Keplerian orbit is projected onto the celestial sphere, the law of 
constant-areal velocity (e.g., Goldstein 1980, Danby 1988, Roy 1988, Murray and Dermott 
1999, Beutler 2004) still holds, where the area is swept by the hne interval between the 
star and the projected common center of mass but not a focus of the observed ellipse. For 
example, let us consider four observed points Pi, P2, P3 and P4 for < ^2 < ^3 < ^4- The 
location {xe,ye) of the projected common center is given (Asada et al. 2004) by 

F1G2 — G1F2 , . 

Xe — —a ^ ^ ^ ^ , (9) 

TP T? T? TP ' ^ ' 

Iliir2 — -ri-C/2 

G1E2-E1G2 , . 

Ve = b-ET-E^ ET-ET, (10) 

where 

= T{j + 1, j) sin Uj+2 + T{j + 2, j + 1) sin Uj 

-T{j + 2,j)smuj+i, (11) 
Fj = T{j + 1, j) cos Uj+2 + T{j + 2, j + 1) cos Uj 

-T{j + 2,j)cosuj+i, (12) 
Gj = T(j + l,j>,+2 + r(j + 2,j + l)M,-T(j + 2,j)w,+i, (13) 
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and the eccentric anomaly in the observed eUipse is given by uj = arctan(a?/j/6xj). What 
Eqs. dn]) and flTU]) tell us is the relative position of the projected common center of mass 
with respect to the center of the observed ellipse, that is the origin of our coordinates 
system. If one wishes to take the origin of the coordinates as another point such as the 
location of the primary star, (xe, ye) must be shifted by the displacement from the center of 
the observed ellipse to the new point. The above solution is generalized to determination of 
an open orbit (Asada 2006). 



In reality, however, we have observational errors in each position measurement. We 
assume that these random errors obey the Gaussian distribution. When the errors vanish, 
any set of four points Pp, Pg, Pf and Pg {tp < tg < < tg) among observational data must 
satisfy Eqs. and ffTU]) . if one replaces 1— >p, 3— >r and 4 s. All of a, b, Ej, Fj 
and Gj in Eqs. ([9]) and (fTOj) have been determined up to this point and Xg and ye are the 
parameters to be estimated. We should note that Eqs. (Q and (ITOl) are linear in parameters 
Xe and ye to be determined by the least square method. Hence, becomes square in Xe 
and ye so that the extremum of can tell us 



where the summation is taken for every set of four points Pp, Pq, Pr and Ps with a 
correspondence a.sj—*p,j + l—*q,j + 2—*r and j + 3 s for defining Ej, Fj and Gj by 
Eqs. (ITTi) - (fT^ . and the number of all the combinations is = n!/4!(n — 4)!. In practice, 
the number of summing, namely a denominator of Eqs. f|T^ and f|T5l) . can be reduced 
significantly from to n/4 in the case of n = Am, if one considers only a set of four points 
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Extension to observational data 




(14) 



(15) 



-8- 



as {Pj, Pj+m, Pj+2m, Pj+3m) for j — 1,2, ■ ■ ■ ,171. It is Worthwhile to mention that each point 
Pj appears only once in this case, and the reduction is useful when applied to a lot of data. 

The mapping between the observed ellipse and the Keplerian orbit parameterized by 
the major axis gk and the ellipticity ex is specified by the inclination angle i and the 
ascending node lu (See Fig. 1). 

Given a, b, Xg and ye, one can analytically determine the remaining parameters e^, i, 
gk and lu in order as (Asada et al. 2004) 



^ + |. (16) 
cosi=iK-v^e-4), (17) 



where 



V 1 + cos^ I 

cos 2a; = — t-^t, (19) 



c = — V^fT^, (20) 



ex 



1 a^yl + b^ 



(c^ + D^)yr^ 

a6 



(22) 



4. Numerical test 

In this section, we verify the present formula by numerical computations for various 
values of the orbital parameters. In order to make a verification, we assume random 

errors in measurements of positions of the star and show that the orbital parameters 
are reconstructed: First, we assume a Keplerian orbit specified by qk and ck and a 
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projection defined by i and u. We prepare twenty seven sets of the parameters as ai^- = 1, 
ex = 0.1, 0.3, 0.6, i = 0, 30, 60 deg. and u = 0, 30, 60 deg. 

For each test, we pick up twelve points in the original Kepler orbit by assuming the 
same time interval between neighboring points for simplicity. Next, for assumed i and u, we 
estimate the projected locations of the twelve points in the apparent ellipse plane. We add 
random errors to each point in order to imitate a "measured position" with observational 
errors (See Fig. 2). In order to produce such a simulated data, we assume that the absolute 
standard deviation of the random error in the position measurement is 0.001 in the units 
of = 1- We apply our formula to the simulated data in order to determine the orbital 
elements, and make a comparison between the assumed value and the retrieved one. Table 
1 shows a good agreement between them. 

For each parameter set, we perform a hundred of numerical runs for a statistical treat- 
ment. We define an error as the absolute standard deviation as Aa = a/ < {a'j^ — axY >, 
Ae = -^/^'(e^^^'exp^, Ai = (i' — iy > and Aa; = < (a;' — a;)^ >, where the prime 
denotes values retrieved by using the present formula and the square bracket denotes a 
mean over each set. In Table 1, asterisk (*) indicates that the formula does not always give 
a real value but a complex one. The case of complex numbers occurs because a large error 
in measurements of the positions in the apparent ellipse plane causes an anomaly in the 
apparent motion such as an apparent clockwise motion even if a true one is anti-clockwise. 

It is worthwhile to make a comment on Table 1: In the case of no inclination as i — 0, 
the ascending node does not exist, so that an error in retrieving cu can be apparently quite 
large, though it is harmless. 
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5. Conclusion 

In this paper, as a simplified method to determine a binary orbit, we extend the 
solution considering observational data. We expect that the present formula will be used for 
orbit determinations by some planned astrometric space missions such as SIMp (Shao 2004), 
GA1A0 (Mignard 2004, Ferryman 2004) and JASMINeS (Gouda et al. 2004), which will 
observe with the accuracy of possibly a few micro arcseconds a number of binaries whose 
components are a visible star and a companion ranging from a black hole to an extra-solar 
planet. Our numerical computations correspond to a case that by these missions we observe 
a binary system at a distance of 1 kpc from us, where the semimajor axis of the primary's 
orbit is around lAU. 

We would like to thank Professor N. Gouda for useful information on future astrometry 
missions. We would like to thank Professor M. Kasai and Professor S. Kuramata for 
continuous encouragement. 



^http: / / sim.jpl.nasa.gov/ 

^http: / / astro.estec.esa.nl/ GAIA/ 

^http: / /www .jasmine-galaxy.org/ 



REFERENCES 

Aitken R. G., 1964 The Binary Stars (NY: Dover) 
Asada H., Akasaka T., Kasai M., 2004, PASJ., 56, L35 



Asada H., 2006, Celest. Mech. m press ( |astro-ph/0609769t ) 
Benedict G. F. et al., 2002, ApJ., 581, LI 15 

Beutler G., 2004 Methods of Celestial Mechanics (Berlin: Springer) 

Catovic Z., Olevic D., 1992 in lAU Colloquim 135, ASP Conference Series Vol. 32 (eds 
McAlister H. A., Hartkopf W. I., ) 217-219 (San Francisco, Astronomical Society of 
the Pacific). 

Danby J. M. A., 1988 Fundamentals of Celestial Mechanics (VA: William-Bell) 

Eichhorn H. K., Xu Y., 1990, ApJ., 358, 575 

Goldstein H., 1980 Classical Mechanics (MA: Addison- Wesley) 

Gouda N. et al., 'Japan Astrometry Satellite Mission for Infrared Exploration (JASMINE)', 
Proc. The Three- Dimensional Universe with Gaia, 4-7 October 2004, Paris 
(Netherlands: ESA Publications) 

Mignard F., 'Overall Science Goals of the Gaia Mission', Proc. The Three- Dimensional 
Universe with Gaia, 4-7 October 2004, Paris (Netherlands: ESA Publications) 

Miyoshi M. et al, 1995, Nature, 373, 127 

Murray C. D., Dermott S. F., 1999 Solar System Dynamics (Cambridge: Cambridge Univ. 
Press) 

Olevic D., Cvetkovic Z., 2004, A&A, 415, 259 



- 12 - 



Ferryman M. A. C, 'Overview of the Gaia Mission', Proc. The Three- Dimensional Universe 
with Gaia, 4-7 October 2004, Paris (Netherlands: ESA Pubhcations) 

Roy A. E., 1988 Orbital Motion (Bristol: Institute of Physics Pubhshing) 

Shao M., 2004, 'Science Overview and Status of the SIM Project', SPIE 5491-36 

Weaver H. A. et al, 2006, Nature, 439, 943 



This manuscript was prepared with the AAS macros v5.2. 



Fig. 1. — A relation between the Keplerian ellipse and the observed one. The thick solid curve 
shows the Keplerian ellipse and the dash-dotted curve shows the observed one projected onto 
a celestial sphere denoted by x-y plane. The line of sight is along the 2;-axis. The common 
center of mass is projected into the point denoted by Pg in the apparent ellipse plane. The 
locations of the intersection of the observed ellipse and the projected major and minor axes 
are denoted by Pl and P5, respectively. 
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Measured Position 




Fig. 2. — A schematic figure of a simulation. The sign "+" denotes a measured position in 
the simulations. The dashed curve is an apparent ellipse. 
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Table 1. Reconstructing the parameters for a simulated data 



o-k ~ ~ i ~ <^ ^1 ^6 [deg.] Au! [deg.] 



1-0.1- 0- 


0.00113 


0.00271 


3.27 


50.9 


1-0.1- 0-30 


0.00106 


0.00258 


3.16 


24.7 


1-0.1- 0-60 


* 


* 


* 


* 


1-0.1-30- 


0.000866 


0.00305 


0.116 


1.34 


1-0.1-30-30 


0.00110 


0.00296 


0.154 


1.30 


1-0.1-30-60 


0.00110 


0.00343 


0.122 


1.04 


1-0.1-60- 


0.00112 


0.00450 


0.0574 


2.08 


1-0.1-60-30 


0.00193 


0.00544 


0.0951 


2.20 


1-0.1-60-60 


0.00141 


0.00501 


0.0605 


1.90 


1-0.3- 0- 


0.00180 


0.00497 


4.29 


47.5 


1-0.3- 0-30 


0.00166 


0.00516 


4.29 


27.8 


1-0.3- 0-60 


0.00193 


0.00555 


4.52 


28.4 


1-0.3-30- 


0.000933 


0.00518 


0.224 


0.943 


1-0.3-30-30 


0.00175 


0.00542 


0.317 


0.719 


1-0.3-30-60 


0.00142 


0.00597 


0.164 


0.449 


1-0.3-60- 


0.00157 


0.00884 


0.122 


1.17 


1-0.3-60-30 


0.00238 


0.00856 


0.150 


0.832 


1-0.3-60-60 


0.00227 


0.00797 


0.0888 


0.715 


1-0.6- 0- 


0.0105 


0.0137 


9.16 


54.6 


1-0.6- 0-30 


0.00977 


0.0147 


9.24 


34.8 


1-0.6- 0-60 


0.0131 


0.0150 


9.48 


33.3 


1-0.6-30- 


0.00240 


0.0168 


1.67 


2.37 


1-0.6-30-30 


0.00374 


0.0172 


1.32 


2.48 


1-0.6-30-60 


0.00953 


0.0150 


0.623 


2.54 


1-0.6-60- 


0.00400 


0.0279 


0.919 


1.68 


1-0.6-60-30 


0.00614 


0.0287 


0.765 


0.966 


1-0.6-60-60 


0.0117 


0.0191 


0.256 


0.586 



